Two-dimensional finite element method
\( \cup_{K \in T_{hp}} = \Omega \), \( \textrm{ meas} K_i \cap K_j = K_i \textrm{ dla } i=j; 0 \textrm{ for } i \neq j \).
\( \forall K \in T_{\frac{h}{2}p+1} \exists K_1,K_2,K_3,K_4 \in {\cal P}(T_{hp},K) \) such that
\( K=K_1\cup K_2 \cup K_3 \cup K_4, \textrm{meas}K_i \cap K_j=K_i \textrm{ for } i=j; 0 \textrm{ for } i \neq j \),
\( X(K_1),X(K_2),X(K_3),X(K_4) \in {\cal P}(T_{hp},X(K)), dim X(K)=\\ dim(K_1)+7=dim(K_2)+7=dim(K_3)+7=dim(K_4)+7 \) where
\( {\cal P}(T_{hp},K) \) and \( {\cal P}(T_{hp},X(K)) \) mean projections on the first and second components \( \left(K, X\left(K\right), \Pi_p \right) \).
The dimension of the space of shape functions over the new four finite elements can be estimated by adding of seven new shape functions, corresponding to increasing the degree of polynomials by 1 inside the element horizontally and vertically (making a total of three additional polynomials - one with the degree increased vertically, the other with the degree increased in the horizontal direction, and the third with the degree increased in both directions), and adding four new shape functions obtained by adding new polynomials on the edges of the elements, with the degree of the polynomial increased by 1.
Approximation space over a two-dimensional computational mesh \( T_{hp}=\{ \left(K, X\left(K\right), \Pi_p \right) \}_K \) is
\( V_{hp }= \{ v \in C(\Omega): \forall K \in {\cal P }( T_{hp},K): {\cal P }( v,K) \in X(K) \} \)
where \( {\cal P }( T_{hp},K) \) is a set of intervals representing the geometry of elements extracted from the triple representing a two-dimensional computational mesh, \( {\cal P }( v,K) \) is a projection of a function onto an interval representing the element's geometry.
Let \( T_{hp}=\{ \left(K, X\left(K\right), \Pi_p \right) \}_K \) mean a two-dimensional computational mesh.
Let \( \{ e_i^{hp} \}_i \) mean the base of space \( V_{hp }= span\{ e_i^{hp} \} \).
Let \( \chi^K_k \in X \left( K\right) \) denote the shape function over the element \( K \).
Then \( \forall K \in {\cal P }( T_{hp},K), \forall i, \exists k : {\cal P}(e_i^{hp},K) = \chi^K_k \).
There is an inverse mapping \( {\cal I}^2 \ni (k,K)\rightarrow i(k,K)\in {\cal I} \) that assigns a number \( i(k,K) \) global i-th base function related to local \( k \)-th shape function over the element \( K \).
In the cases discussed in this manual, this mapping is isomorphic.
Approximation space above a two-dimensional reference mesh \( T_{\frac{h}{2}p+1}=\{ \left(K, X\left(K\right), \Pi_p \right) \}_K \) is
\( V_{\frac{h}{2}p+1 }= \{ v \in C(\Omega): \forall K \in {\cal P }( T_{\frac{h}{2}p+1},K): {\cal P }( v,K) \in X(K) \} \)
where \( {\cal P }( T_{\frac{h}{2}p+1},K) \) is a set of intervals representing the geometry of elements extracted from the triple representing a two-dimensional computational mesh, \( {\cal P }( v,K) \) is a projection of a function onto an interval representing the element's geometry.
Let \( T_{\frac{h}{2}p+1}=\{ \left(K, X\left(K\right), \Pi_p \right) \}_K \) stands for a two-dimensional reference mesh.
Let \( \{ e_i^{\frac{h}{2}p+1} \}_i \) means the base of space \( V_{\frac{h}{2}p+1 }= span\{ e_i^{\frac{h}{2}p+1} \} \).
Let \( \chi^K_k \in X \left( K\right) \) denotes the shape function over the element \( K \).
Then \( \forall K \in {\cal P }( T_{\frac{h}{2}p+1},K), \forall i, \exists k : {\cal P}(e_i^{\frac{h}{2}p+1},K) = \chi^K_k \).
There is an inverse mapping \( {\cal I}^2 \ni (k,K)\rightarrow i(k,K)\in {\cal I} \) which assigns a global index \( i(k,K) \) of base function based on the local index of \( k \)-th shape function over the element \( K \).
In the cases discussed in this manual, this mapping is isomorphic.
Approximation space above a two-dimensional reference computational mesh \( T_{\frac{h}{2}p+1}=\{ \left(K, X\left(K\right), \Pi_p \right) \}_K \) to \( V_{\frac{h}{2}p+1}=span \{ e_j^{\frac{h}{2}p+1} : \forall K \in {\cal P}(T_{\frac{h }{2}p+1},K), \forall \phi_j \in X(K), \exists ! e_i^{\frac{h}{2}p+1} :{\cal P }( e_i^{\frac{h}{2}p+1 },K)=\psi_k \} \) where \( e_i^{\frac{h}{2}p+1} \) is a global base function (element of \( V_{\frac{h}{2}p+1} \) ), and \( \psi_k \) is the local shape function above the element \( K \),
\( {\cal I}^2 \ni (k,K)\rightarrow i(k,K)\in {\cal I} \) is a mapping that assigns a number \( i(k,K) \) global base function index related to local \( k \)-th shape function over the element \( K \).
The reference mesh is used to estimate the relative error over the computational mesh. Often, when a computational mesh is described in the context of a reference mesh, it is called a sparse mesh (or coarse mesh), and a reference mesh is called a fine mesh (or a dense mesh), because it is created by breaking the elements into smaller elements and increasing the polynomial order of approximation by one.
For a given computational mesh \( T_{hp}=\{ \left(K, X\left(K\right), \Pi_p \right) \}_K \) find the coefficients \( \{ u^{hp}_i \}_{i=1,...,N^{hp}} \) of the approximate solution \( V \supset V_{hp } \ni u_{hp }=\sum_{i=1,...,N^{hp } } u_i^{hp } e_i^{hp } \) such that \( \sum_{m=1,...,N^{hp } } u_m^{hp } B(e_m^{hp},e_n^{hp})=L(e_n^{hp}) n=1,...,N^{hp} \) where
\( B(e_m^{hp},e_n^{hp})= \int_{\Omega} \left( \sum_{i=1,2} \sum_{j=1,2} a_{ij}(x_1,x_2) \frac{\partial e_m^{hp}(x_1,x_2) }{\partial x_j } \frac{\partial e_n^{hp}(x_1,x_2) }{\partial x_j }\right. \\ \left. + \sum_{j=1,2} b_j(x_1,x_2) \frac{\partial e_m^{hp}(x_1,x_2) }{\partial x_j } e_n^{hp}(x_1,x_2) +c (x_1,x_2 )e_n^{hp}(x_1,x_2) \right) dx_1 dx_2 + \)
\( +\int_{\Gamma_R} \beta(x_1,x_2) e_m^{hp}(x_1,x_2) e_n^{hp}(x_1,x_2) ds \)
\( L(e_n^{hp})= \int_{\Omega} f(x_1,x_2)e_n^{hp}(x_1,x_2) dx_1 dx_2 + \int_{\Gamma_R} g(x_1,x_2) e_n^{hp}(x_1,x_2) ds \).
Base \( \{ e^{hp}_i \}_{i=1,...,N^{hp}} \) of the approximation space \( V_{hp} \) is obtained by summing shape functions from space \( X\left(K\right) = {\cal P}(T_{hp},X(K)) \) for particular elements from the mesh \( T_{hp} \) into global base functions.
For a given computational mesh \( T_{\frac{h}{2}p+1}=\{ \left(K, X\left(K\right), \Pi_p \right) \}_K \) find the coefficients \( \{ u^{\frac{h}{2}p+1}_i \}_{i=1,...,n^{\frac{h}{2}p+1}} \) of approximate solution \( V \supset V_{\frac{h}{2}p+1 } \ni u_{\frac{h}{2}p+1 }=\sum_{i=1,...,N^{\frac{h}{2}p+1 } } u_i^{\frac{h}{2}p+1 } e_i^{\frac{h}{2}p+1 } \) such that \( \sum_{m=1,...,N^{\frac{h}{2}p+1 } } u_m^{\frac{h}{2}p+1 } B(e_m^{\frac{h}{2}p+1},e_n^{\frac{h}{2}p+1})=L(e_n^{\frac{h}{2}p+1}) n=1,...,N^{hp} \)
\( B(e_m^{\frac{h}{2}p+1},e_n^{\frac{h}{2}p+1})= \int_{\Omega} \left( \sum_{i=1,2} \sum_{j=1,2} a_{ij}(x_1,x_2) \frac{\partial e_m^{\frac{h}{2}p+1}(x_1,x_2) }{\partial x_j } \frac{\partial e_n^{\frac{h}{2}p+1}(x_1,x_2) }{\partial x_j }\right. \\ \left. + \sum_{j=1,2} b_j(x_1,x_2) \frac{\partial e_m^{\frac{h}{2}p+1}(x_1,x_2) }{\partial x_j } e_n^{\frac{h}{2}p+1}(x_1,x_2) +c (x_1,x_2 )e_n^{\frac{h}{2}p+1}(x_1,x_2) \right) dx_1 dx_2 + \)
\( +\int_{\Gamma_R} \beta(x_1,x_2) e_m^{\frac{h}{2}p+1}(x_1,x_2) e_n^{\frac{h}{2}p+1}(x_1,x_2) ds \)
\( L(e_n^{\frac{h}{2}p+1})= \int_{\Omega} f(x_1,x_2)e_n^{\frac{h}{2}p+1}(x_1,x_2) dx_1 dx_2 + \int_{\Gamma_R} g(x_1,x_2) e_n^{\frac{h}{2}p+1}(x_1,x_2) ds \).
Base \( \{ e^{\frac{h}{2}p+1}_i \}_{i=1,...,N^{\frac{h}{2}p+1}} \) of the approximation space \( V_{hp} \) is obtained by summing shape functions from space \( X\left(K\right) = {\cal P}(T_{\frac{h}{2}p+1},X(K)) \) for individual mesh elements \( T_{\frac{h}{2}p+1} \) into global basis functions.
Coefficients \( \{ u^{hp}_i \}_{i=1,...,N^{hp}} \) is calculated by solving a system of equations \( {\bf Ax}={\bf b} \) where
\( {\bf A}= \begin{bmatrix} B(e_1^{hp},e_1^{hp}) & \cdots & B(e_{N^{hp}}^{hp},e_1^{hp}) \\ \vdots & \cdots & \vdots \\ B(e_{N^{hp}}^{hp},e_1^{hp}) & \cdots & B(e_{N^{hp}}^{hp},e_{N^{hp}}^{hp}) \end{bmatrix} \) ,
\( {\bf x}= \begin{bmatrix} u_1^{hp} \\ \vdots\\ u_{N^{hp}}^{hp} \end{bmatrix} \),
\( {\bf b}= \begin{bmatrix} L(e_1^{hp}) \\ \vdots \\ L(e_{N^{hp}}^{hp}) \end{bmatrix} \).
The Dirichlet boundary condition is obtained by setting to zero the rows corresponding to the basis functions defined on the Dirichlet boundary, placing one on the diagonal and setting zero on the right side. The effect of the boundary condition on our system of equations has been transferred to the right-hand side by a shift operation.
Integrals \( B(e_m^{hp},e_n^{hp}) \) and \( L(e_n^{hp}) \) are computed element by element using the appropriate Gaussian quadratures.
\( B(e_m^{hp},e_n^{hp}) = \sum_k \left( \sum_{i=1,2} \sum_{j=1,2} a_{ij}(x^k_1,x^k_2) \frac{\partial e_m^{hp}(x^k_1,x^k_2) }{\partial x_j } \frac{\partial e_n^{hp}(x^k_1,x^k_2) }{\partial x_j } \right. \\ \left. + \sum_{j=1,2} b_j(x^k_1,x^k_2) \frac{\partial e_m^{hp}(x^k_1,x^k_2) }{\partial x_j } e_n^{hp}(x^k_1,x^k_2) +c (x^k_1,x^k_2 )e_n^{hp}(x^k_1,x^k_2) \right) Jac(x^k_1,x^k_2) w^k_1 w^k_2 + \)
\( + \sum_k \beta(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) e_m^{hp}(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) e_n^{hp}(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) Jac(z^k) w^k \)
\( L(e_n^{hp})= \sum_k f(x^k_1,x^k_2)e_n^{hp}(x^k_1,x^k_2) Jac(x^k_1,x^k_2)w^k + \sum_k g(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) e_n^{hp}(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) Jac(z^k)w^k \) where
\( [0,1] \ni z \rightarrow (x^{\Gamma_R}_1(z),x^{\Gamma_R}_2(z))\in \Gamma_R \) is a function that parameterizes the Robin boundary.
Generation of a system of equations for a given computational mesh \( T_{hp}=\{ \left(K, X\left(K\right), \Pi_p \right) \}_K \)
can be obtained with the following algorithm:
1 for \( K \in {\cal P}(T_{hp },K) \) (loop through elements of the computational mesh)
2 for \( \psi_k \in X(K) \) (loop through shape functions on an element \( K \) )
3 \( {\bf b}(i(k,K)) += L(\psi_k) \) (computing the local integral contribution to an element \( K \) )
4 for \( \psi_l \in X(K) \) (loop through shape functions on an element \( K \) )
5 \( {\bf A}(i(k,K),i(l,K)) += B(\psi_k,\psi_l) \) (computing the local integral contribution to an element \( K \) ).
Annotation 2: Solving the problem of the finite element method over a two-dimensional reference mesh
Coefficients \( \{ u^{\frac{h}{2}p+1}_i \}_{i=1,...,N^{hp}} \) is calculated by solving a system of equations \( {\bf Ax}={\bf b} \) where \( {\bf A}= \begin{bmatrix} B(e_1^{\frac{h}{2}p+1},e_1^{\frac{h}{2}p+1}) & \cdots & B(e_{N^{\frac{h}{2}p+1}}^{\frac{h}{2}p+1},e_1^{\frac{h}{2}p+1}) \\ \vdots & \cdots & \vdots \\ B(e_1^{\frac{h}{2}p+1},e_{N^{\frac{h}{2}p+1}}^{\frac{h}{2}p+1}) & \cdots & B(e_{N^{\frac{h}{2}p+1}}^{\frac{h}{2}p+1},e_{N^{\frac{h}{2}p+1}}^{\frac{h}{2}p+1}) \end{bmatrix} \) ,
\( {\bf x}= \begin{bmatrix} u_1^{\frac{h}{2}p+1} \\ \vdots\\ u_{N^{\frac{h}{2}p+1}}^{\frac{h}{2}p+1} \end{bmatrix} \),
\( {\bf b}= \begin{bmatrix} L(e_1^{\frac{h}{2}p+1}) \\ \vdots \\ L(e_{N^{\frac{h}{2}p+1}}^{\frac{h}{2}p+1}) \end{bmatrix} \).
The Dirichlet boundary condition is obtained by setting to zero the rows corresponding to the basis functions defined on the Dirichlet boundary, placing one on the diagonal and zero in the right-hand side. The effect of the boundary condition on our system of equations has been transferred to the right-hand side by the shift operation.
Integrals \( B(e_m^{\frac{h}{2}p+1},e_n^{\frac{h}{2}p+1}) \) and \( L(e_n^{hp}) \) are computed using the appropriate Gaussian quadratures.
\( B(e_m^{\frac{h}{2}p+1},e_n^{\frac{h}{2}p+1}) = \sum_k \left( \sum_{i=1,2} \sum_{j=1,2} a_{ij}(x^k_1,x^k_2) \frac{\partial e_m^{\frac{h}{2}p+1}(x^k_1,x^k_2) }{\partial x_j } \frac{\partial e_n^{\frac{h}{2}p+1}(x^k_1,x^k_2) }{\partial x_j } \right. \\ \left. + \sum_{j=1,2} b_j(x^k_1,x^k_2) \frac{\partial e_m^{\frac{h}{2}p+1}(x^k_1,x^k_2) }{\partial x_j } e_n^{\frac{h}{2}p+1}(x^k_1,x^k_2) +c (x^k_1,x^k_2 )e_n^{\frac{h}{2}p+1}(x^k_1,x^k_2) \right) Jac(x^k_1,x^k_2) w^k_1 w^k_2 + \)
\( + \sum_k \beta(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) e_m^{hp}(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) e_n^{\frac{h}{2}p+1}(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) Jac(z^k) w^k \)
\( L(e_n^{hp})= \sum_k f(x^k_1,x^k_2)e_n^{\frac{h}{2}p+1}(x^k_1,x^k_2) Jac(x^k_1,x^k_2)w^k + \sum_k g(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) e_n^{\frac{h}{2}p+1}(x^{\Gamma_R}_1(z^k),x^{\Gamma_R}_2(z^k)) Jac(z^k)w^k \) where
\( Jac(x^k) \) is the value at the Gaussian Jacobian quadrature point of the mapping from the reference element to the given element, \( [0,1] \ni z \rightarrow (x^{\Gamma_R}_1(z),x^{\Gamma_R}_2(z))\in \Gamma_R \) is a function that parameterizes the Robin boundary.
Generation of a system of equations for a given calculation mesh \( T_{\frac{h}{2}p+1}=\{ \left(K, X\left(K\right), \Pi_p \right) \}_K \)
can be obtained with the following algorithm:
1 for \( K \in {\cal P}(T_{\frac{h}{2}p+1 },K) \) (loop through computational mesh elements)
2 for \( \psi_k \in X(K) \) (loop through shape functions on an element \( K \) )
3 \( {\bf b}(i(k,K)) += L(\psi_k) \) (computing the local integral contribution to an element \( K \) )
4 for \( \psi_l \in X(K) \) (loop through shape functions on an element \( K \) )
5 \( {\bf A}(i(k,K),i(l,K)) += B(\psi_k,\psi_l) \) (computing the local integral contribution to an element \( K \) ).